In [1]:
%pylab inline
import numpy as np
import numpy.random as rand
from simpledist import distributions as dists #https://github.com/timothydmorton/simpledist
from starutils.populations import BinaryPopulation, StarPopulation #https://github.com/timothydmorton/starutils
from isochrones import dartmouth #https://github.com/timothydmorton/isochrones
dar = dartmouth.Dartmouth_Isochrone()
In [2]:
lo_mass = 0.15
hi_mass = 0.55
imf = dists.PowerLaw_Distribution(-1.3,lo_mass,hi_mass) #kroupa IMF with alpha=-1.3 in this mass range
Simulate primary and secondary stellar populations:
In [3]:
n = 1e5
feh = 0.0
age = 9.6
#figure out masses
minmass = 0.11
m1s = imf.rvs(n)
qmin = minmass/m1s
qs = rand.random(n)*(1-qmin) + qmin #flat mass ratio distribution
m2s = m1s*qs
#stellar properties from isochrones
primaries = dar(m1s,age,feh) #generates populations from Dartmouth models
secondaries = dar(m2s,age,feh)
#volume-limited sample out to 500pc
max_distance = 500
distance_distribution = dists.PowerLaw_Distribution(2.,1.,max_distance) #p(d) ~ d^2
distances = distance_distribution.rvs(n)
binaries = BinaryPopulation(primaries,secondaries,distance=distances)
singles = StarPopulation(primaries,distance=distances)
In [4]:
binaries.prophist('Kepler_mag',histtype='step',label='binaries')
singles.prophist('Kepler_mag',histtype='step',fig=0,label='singles')
legend(loc='upper left')
axvline(16,color='r',ls=':')
title('Volume limited sample out to 500pc')
Out[4]:
In [5]:
limiting_mag = 16
binaries.constrain_property('Kepler_mag',hi=limiting_mag)
singles.constrain_property('Kepler_mag',hi=limiting_mag)
In [6]:
binaries.constraint_piechart(title='Binaries')
singles.constraint_piechart(title='Singles')
We started with equal numbers of binaries and singles, i.e. total binary fraction of 50%. Now, the binary fraction in the selected sample is:
In [7]:
binaries.selectfrac/(binaries.selectfrac + singles.selectfrac)
Out[7]:
Now, if the original binary fraction was only 30%, with 70% in single systems, the new binary fraction would be (I think this is the correct calculation):
In [8]:
(0.3/0.5)*binaries.selectfrac/((0.3/0.5)*binaries.selectfrac + (0.7/0.5)*singles.selectfrac)
Out[8]:
--a 1.4x increase of the binary fraction.
Now, to check out what happens with transit probability:
In [9]:
selected_rsum = np.array(binaries.selected['radius_A'] + binaries.selected['radius_B'])
original_rsum = np.array(binaries.stars['radius_A'] + binaries.stars['radius_B'])
selected_msum = np.array(binaries.selected['radius_A'] + binaries.selected['radius_B'])
original_msum = np.array(binaries.stars['radius_A'] + binaries.stars['radius_B'])
selected = selected_rsum/selected_msum**(1./3)
original = original_rsum/original_msum**(1./3)
hist(selected,histtype='step', normed=True, label='magnitude selected')
hist(original,histtype='step', normed=True, label='volume limited')
legend()
axvline(selected.mean(),color='b',ls=':')
axvline(original.mean(),color='g',ls=':')
xlabel('relative transit probability $(R_1+R_2)/(M_1 + M_2)^{1/3}$');
Bottom line: An intial binary fraction of 50% becomes 62% in the magnitude limited sample. And this selection also favors larger transit probabilities by a factor of 25-30%.
Another issue that hasn't yet been quantified here: very close binaries are very often found in triple systems, which will further exacerbate these selection effects; that is, triple systems will be even more significantly over-represented compared to their intrinsic occurrence rates, given the magnitude selection, and nearly all of these will be hierarchical, with two having short periods.
OK, trying this same exercise with new class defined that packages up the set-up calculations above:
In [8]:
n = 1e5
feh = rand.normal(1e5)*0.1 - 0.1
age = 9.6
from starutils.populations import VolumeLimitedPopulation
m1s = imf.rvs(n)
m1s = np.ones(n)*0.5
pop2 = VolumeLimitedPopulation(m1s, 500, binary_fraction=0.3, n=n)
In [9]:
print pop2.binary_fraction()
print pop2.binary_fraction('Kepler_mag < 16')
In [4]:
from starutils.populations import MultipleStarPopulation
In [16]:
f2 = 0.3; f3 = 0.08
s2 = 0.1; s3 = 0.7
pop = MultipleStarPopulation(.5, max_distance=500, f_binary=f2, f_triple=f3)
In [17]:
b2 = pop.binary_fraction('Kepler_mag < 16')/pop.binary_fraction()
b3 = pop.triple_fraction('Kepler_mag < 16')/pop.triple_fraction()
In [18]:
print s2*f2 + s3*f3
print s2*f2*b2 + s3*f3*b3
In [ ]: